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ABSTRACT 

We present the star formation histories, luminosities, colors, mass to light ratios, and halo 
masses of "galaxies" formed in a simulation of cosmological reionization. We compare these 
galaxies with Lyman Break Galaxies observed at high redshift. While the simulation is severely 
limited by the small box size, the simulated galaxies do appear as fainter cousins of the observed 
ones. They have similar colors, and both simulated and observed galaxies can be fitted with the 
same luminosity function. This is significant to the process of reionization since in the simulation 
the stars alone are responsible for reionization at a redshift consistent with that of the Gunn- 
Peterson trough seen in the absorption spectra of QSO's. The brightest galaxies at z — 4 are 
indeed quite young in accord with observational studies. But these brightest galaxies are free 
riders - they contributed only about 25% to the reionization of the universe at z > 6. Instead, the 
bulk of the work was done by dimmer galaxies, those that fall within the 28 < I < 30 magnitude 
range at z ~ 4. These dimmer galaxies are not necessarily less massive than the brightest ones. 

Subject headings: cosmology: theory - cosmology: large-scale structure of universe - galaxies: evolution 
- galaxies: formation - galaxies: high-redshift - galaxies: starburst - galaxies: stellar content - infrared: 
galaxies 



1. Introduction 

The most dramatic global transition in the re- 
cent history of the universe is the reionization of 
the intergalactic medium. The recent identifica- 
tion of the predicted Gunn-Peterson trough in the 
absorption spectra of several QSO's has placed 
this event at about a redshift of 6 (Djorgovski et 
al. 2002; Becker et al. 2001). 

Ionizing radiation from newly formed stars has 
been suggested as the most likely source of this 
radiation. Bright, optically selected QSO's ap- 
pear to be too sparse during the relevant period 
(Fan et al. 2001). Low brightness AGN's remain 
a possibility, but no other strong candidates have 
emerged. 

It is natural to look for evidence for these ion- 
izing stars in the accumulating observational data 
on high redshift (z < 4) galaxies (Madau et al. 
1996; Madau, Pozzetti, & Dickinson 1998; Stei- 
del et al. 1999; Ouchi et al. 2001; Ferguson, Dick- 
inson, k Papovich 2002; Hu et al. 2002; Malho- 



tra k Rhoads 2002; Lehnert k Bremer 2003; Ko- 
daira et al. 2003; Iwata et al. 2003). Some studies 
on Lyman Break Galaxies (LBG) have questioned 
whether enough stars could have been formed to 
accomplish reionization at so early a time (Fergu- 
son, Dickinson, k Papovich 2002; Dickinson et al. 
2003). 

In the popular Lambda-CDM cosmological 
model the seeds of galaxy formation lie in minute 
density fluctuations, which manifest themselves 
as temperature fluctuations in the Cosmic Mi- 
crowave Background (CMB). As the fluctuations 
continue to grow into dark matter halos, they trap 
gas within their potential wells. That part of the 
gas that is able to radiate its energy efficiently 
eventually condenses into dense molecular clouds, 
that give birth to stars. 

Reionization has been studied within this 
paradigm both by semi-analytical methods (Giroux 
k Shapiro 1996; Ciardi k Ferrara 1997; Haiman 
k Loeb 1997, 1998; Madau et al. 1999; Valageas k 
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Silk 1999; Chili & Ostriker 2000; Miralda-Escude, 
Haehnelt & Rees 2000; Barkana & Loeb 2000; 
Oh ct al. 2001; Wyithe & Loeb 2003) and by 
numerical simulations (Ostriker & Gnedin 1997; 
Gnedin & Ostriker 1998; Gnedin 2000; Nakamoto, 
Umemura, & Susa 2001; Razoumov et al. 2002). 
The former necessarily involve broad simplifying 
assumptions and, therefore, have limited predic- 
tive power. The latter, although more realistic and 
fully self-consistent, are limited by the dynamical 
range that can be accommodated by existing nu- 
merical techniques on modern supercomputers. 
The two approaches are complementary and will 
have to remain so until the next generation of 
cosmological numerical codes become available on 
faster computers. 

For this paper we use a simulation that was re- 
ported in Gnedin (2000). The simulation includes 
dark matter, gas, star formation, chemistry and 
ionization balance in the primordial plasma, and 
an approximate implementation of 3D radiative 
transfer. All galaxies with total masses in excess 
of 3 x 1O 8 M are resolved both spatially and by 
mass. The simulation was continued until z = 4.0. 

The overall picture for the evolution of the early 
universe that has emerged from this simulation is 
as follows: 

Star formation first begins in earnest at about 
200 Myr (z ~ 20). The overall star formation rate 
increases to a maximum at about 1 Gyr (z ~ 5.5) 
and decreases slightly thereafter. 

Rcionization of the intergalactic medium oc- 
curs over a prolonged period of time from about 
350 Myr (z ~ 12) to about 900 Myr (z ~ 6), with 
expanding Hn regions in the low density IGM 
merging at z ~ 7. This latter moment of rapid 
transition is often identified as the "moment of 
reionization" . This moment approximately cor- 
responds to the mass (volume) averaged neutral 
fraction in the low density IGM falling below 1 
(0.1) percent respectively. 

As a result of rcionization, the gas content of 
halos less massive than about 3 x 10 9 M Q solar 
masses is drastically reduced. 

Since this simulation does not include AGN's, 
it provides a detailed example of how stars alone 
might be sufficient to accomplish reionization. 
The addition of AGN's to this scenario would 
presumably lead to even more ionizing photons, 



which could take up the slack in the event that the 
amount of star formation in the simulation were 
overestimated. 

The simulation provides us with a detailed star 
formation history for each galaxy and thus allows 
us to make detailed photometric predictions. In 
this paper we present the star formation histories, 
luminosities, colors, mass to light ratios, and halo 
masses of the simulated galaxies. 

Our simulation, while state of the art, is only 
barely sufficient for our purpose, because the com- 
putational box is so small. This inherent limita- 
tion of any simulation limits the number of bright 
objects that can be modeled, while the obser- 
vations approach the galaxy hierarchy from the 
other end, detecting bright galaxies and missing 
the dim ones. In this sense we can think of ob- 
servations and simulations as moving toward each 
other along the galaxy luminosity curve. Although 
there is not much overlap yet, the junction point 
has at last been reached, as we show below. 

We argue that to the extent comparisons can 
be made, our predictions are consistent with ob- 
servations of Lyman Break Galaxies. The varied 
star formation histories that we see provide the key 
to understanding why the observations appear to 
argue that LBG galaxies are too young to have 
reionized the universe. 

2. Method 
2.1. Simulation 

The specific simulation we use in this paper 
is fully described in Gnedin (2000). The simula- 
tion was performed using the Softened Lagrangian 
Hydrodynamics (SLH) code, which includes dark 
matter, gas, star formation, chemistry and ion- 
ization balance in the primordial plasma, and 3D 
radiative transfer (in an approximate implemen- 
tation). All these physical ingredients are re- 
quired to properly model the process of cosmo- 
logical reionization. 

The simulation of a representative CDM+A 
cosmological model 1 was performed in a comov- 
ing box with the size of 4/i _1 Mpc with the total 

With the following cosmological parameters: f!o = 0.3, 
Q A = 0.7, h = 0.7, Q b = 0.04, n = 1, <r 8 = 0.91, where 
the amplitude and the slope of the primordial spectrum are 
fixed by the COBE and cluster-scale normalizations. 
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mass resolution of 4 x 10 6 M and the comoving 
spatial resolution of l/i -1 kpc. This resolution is a 
reasonable compromise between the need to have 
the box large enough to accommodate several Hn 
regions and the need to have spatial resolution 
high enough to resolve sources of ionizing radi- 
ation. Convergence studies presented in Gncdin 
(2000) demonstrate that, while still limited, the 
resolution of this simulation is sufficient to model 
reionization at a semi-qualitative level (with pre- 
cision of the order of 50% or so). 

The simulation was stopped at z = 4 because at 
this time the rms density fluctuation in the com- 
putational box is about 0.25, and at later times 
the box ceases to be a representative region of the 
universe even for the dark matter. 

The simulation was designed in such a way 
as to approximately reproduce the measured 
star formation rate density at z = 4 of about 
0.1M Q /yr/Mpc 3 (Steidel et al. 1999). 

Our simulation assumes that all cosmologi- 
cal reionization occurs via radiation from stel- 
lar sources, with star formation parameterized by 
the phenomenological Schmidt law as discussed in 
Gnedin (2000a). An alternative and complemen- 
tary scenario would be one in which the bulk of 
reionization is produced by Active Galactic Nu- 
clei (AGN) . In light of recent high-redshift quasar 
counts (e.g., Fan et al. 2001), a scenario in which 
the universe is reionized by bright optically se- 
lected QSOs seems unlikely in any event. There 
remains however the possibility that low bright- 
ness AGNs contributed significantly (if not domi- 
nantly) to the reionization of the universe. 

For the purpose of this paper, however, this 
simulation is quite suitable for the following rea- 
son. Our main goal in this paper is to investi- 
gate whether galaxies are capable of reionizing the 
universe and still appear as young as they do at 
z ~ 4. Thus, by ignoring the AGN component, 
we consider the most extreme case, since we re- 
quire galaxies to carry all the burden of reioniz- 
ing the universe by themselves. One would then 
expect that if galaxies are capable of reionizing 
the universe and appear as young as they do at 
z ~ 4, then it will be even easier for the models to 
agree with observations if some ionization is done 
by AGNs. 

We identify galaxies in the simulation with 



gravitationally bound objects, which are selected 
with the DENMAX algorithm of Bertschinger & 
Gelb (1991). 

2.2. Population Synthesis 

Population synthesis is carried out using the 
Starburst99 package (Leitherer et al. 1999). The 
version we use included the July 2002 update, 
which revised the UV spectra for O and B stars. 
These spectra, based on the recent work of Smith, 
Norris, & Crawther (2002), have decreased ioniz- 
ing radiation during the early stages of star forma- 
tion. These changes affect not only the amount of 
Lyman alpha continuum, but also the amount of 
Lyman alpha that can be produced if this UV is 
absorbed by local gas. We choose the instanta- 
neous burst option. This is appropriate because 
we apply spectral synthesis to individual stellar 
particles formed in the simulation, which sample 
the stellar distribution function at discrete time 
steps. 

We choose the lowest metallicity option (5% 
solar) available in Starburst99, since the simu- 
lated galaxies have comparable metallicities. Star- 
burst99 maintains same metallicity for each sub- 
sequent generation of stars. This simplification 
should be sufficient for times less than about 1 Gyr 
(Leitherer et al. 1999). 

Ultra-violet absorption by the IGM material 
along the line of sight (i.e. by the Lyman alpha 
forest) is taken into account using the method of 
Madau (1995). 

Reddening corrections are made by the method 
of Calzetti (1997), including the recent calibration 
data of Leitherer et al. (2002) We find it sufficient 
to consider a range of reddening corresponding to 
a color excess of E(B — V) = — 0.3. Higher 
amounts of reddening produce colors inconsistent 
with nearly all the observed galaxies. This range 
is in line with the estimate of Shapely et al. (2001) 
of an average color excess of E{B — V) = 0.15 for 
a sample of Lyman break galaxies at a redshift 
of about 3. This range is also generally consistent 
with the results of others for high redshift galaxies. 

We also allow for the partial escape of ionizing 
photons from the modeled galaxies by reducing 
the flux above 13.6 eV by a factor /esc, an d con- 
verting 2/3 of the removed photons into Lyman- 
alpha emission following the recipe of Malhotra & 
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Fig. 1. — Calculated spectrum of a typical bright 
simulated galaxy at redshift 4.0. Ultra-violet ab- 
sorption by the Lyman alpha forest has been taken 
into account. An escape fraction of 1.0 was used 
with no reddening correction. 

Rhoads (2002). The escape fraction /esc becomes 
another parameter in our modeling. We consider 
three values for this parameter: 0.1, 0.5, and 1.0. 

The AB magnitude system is used through- 
out. We adopt the same cosmology as used in 
the simulation (flat concordance cosmology with 
H = 70km/s/Mpc. 

For comparison purposes the G — I colors for 
observed galaxies over a range of rcdshifts were 
corrected to a redshift of 4.0. This was done as 
follows. For each of a range of z values, including 
4.0, a series of 199 model galaxies of different ages 
was constructed at 5 Myr intervals, each having a 
single burst of star formation. For each observed 
galaxy the set of models at the nearest z value was 
selected. The two bracketing models at this red- 
shift plus the two corresponding models at redshift 
4.0 were used to obtain a corrected color using a 
linear interpolation. This method was feasible be- 
cause the relationship between G — I and age is, 
with few exceptions at these intervals, monotonic. 
A similar correction for 1Z — I was not possible. 

Figure 1 shows an illustrative spectrum of a 
typical bright galaxy from our simulation. Sharp 
absorption edges due to intervening Lyman-a and 
Lyman-/? forests are clearly visible. 



Fig. 2. — Luminosity Function at rest frame 1600 
A for simulated galaxies at z = 4.0 in the I band. 
Squares show the unreddened luminosity function, 
while triangles give the luminosity function with 
E(B — V) = 0.15 reddening included. Vertical 
error-bars are 1-sigma Poisson errors. Diamonds 
show the observed luminosity function of HDF 
galaxies from Madau et al. (1996), while crosses 
and asterisks mark the observed LBG luminos- 
ity functions at z ~ 4 and z ~ 3 respectively 
from Steidel et al. (1999). The dashed line is our 
Schechter function fit to the combined luminosity 
function. 

3. Results 

Figure 2 shows a computed luminosity function 
for the simulated galaxies at z = 4.0. As a mea- 
sure of luminosity we use the / band AB mag- 
nitude, which corresponds to a rest frame wave- 
length of about 1600 A. Also shown for compari- 
son is the observed luminosity function for Lyman 
Break Galaxies at z ~ 3 — 4 from a variety of 
authors (Madau et al. 1996; Steidel et al. 1999). 
Because of the difficulties in observing at these 
redshifts, Steidel et al. (1999) have not attempted 
to make a Schechter fit at z = 4, but instead have 
observed that the data are well described by the 
Schechter function fitted to Lyman Break Galax- 
ies at a redshift of about 3. This fit, however, 
does not match the slope of the luminosity func- 
tion of simulated galaxies at faint magnitudes, but 
the combined luminosity function (simulated plus 
observed) can be fit by a single Schechter function 
with a shallower slope. The parameters of both 
fits arc listed in Table 1. 
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Because of the small size of the simulation box, 
only a few of the observed points are expected 
to overlap the simulated ones, and these are all 
from the Hubble Deep Field North, which probed 
fainter magnitudes than did the ground based ob- 
servations of Steidel et al. (1999). But it appears 
that both the simulation and the observations can 
be fitted with the same luminosity function, im- 
plying that the simulation produces galaxies sim- 
ilar to the real ones. The limitation of the small 
box size directly translates into the upper limit on 
galaxy luminosities (there should be at least a few 
galaxies of luminosity L in the simulation box to 
determine 4>{L)). As one can see, our simulation 
has just enough volume to reach the luminosities 
of HDF galaxies, and is a factor of 2 too small 
to have a single galaxy as large as those found 
by Steidel et al. (1999) in the simulation volume. 
Thus, we are not able to make comparisons on a 
galaxy- by-galaxy basis yet, but future simulation 
will be able to model volumes comparable to those 
covered by observations in only a few years. 

Given a Schechter function fit, we can integrate 
the total luminosity function to estimate the frac- 
tion of starlight in observed and simulated galax- 
ies. We find that the observed galaxies contain 
about 2/3 of the total luminosity density, while 
the simulated ones account for about 1/3. The 
contribution of observed galaxies to reionization is 
expected to be somewhat smaller however (Gnedin 
2000b). This is because these galaxies tend to be 
clustered (Steidel et al. 1999; Adelberger et al. 
1998), and most of the time several of them sit 
inside the same Hn region. When an Hn region 
grows to be larger than the mean free path of ion- 
izing photons, these photons are wasted. 

These arguments imply that the rcdshift of 
reionization in the simulation is an underestimate, 
but not by a large factor because the star forma- 
tion rate in the simulation increases rapidly with 
time (Gnedin 2000a). If we assume that the sim- 
ulation underestimates the volume filling fraction 
of the ionized gas at any given time by, say, 50%, 



Table 1: Schechter function fit parameters 



Parameter 


Steidel et al. fit 


Our fit 


a 


-1.60 


-1.22 


M* 


-20.9 


-20.4 


</>* (Mpc~ 3 ) 


1.8 x 10~ 2 


5.7 x 10~ 2 



then the redshift of reionization is underestimated 
in the simulation by about Az ~ 1.5. As we em- 
phasized above, our simulation is clearly insuffi- 
cient to give a quantitatively accurate model of 
reionization, so this level of precision is more than 
acceptable for our purposes. 

A remarkable feature of Fig. 2 is flattening of 
the luminosity function in the dwarf galaxy range. 
This effect is the result of the inhibited accretion 
of gas onto low mass halos because of photoion- 
ization, sometimes improperly called "photoevap- 
oration" (Thoul & Weinberg 1996; Quinn, Katz, 
& Efstathiou 1996; Weinberg, Hernquist, & Katz 
1997; Navarro & Steinmetz 1997; Gnedin 2000b; 
Chiu, Gnedin, & Ostriker 2000; Sommcrville 2002; 
Benson et al. 2002a, b), and has been reproduced 
by earlier simulations (Nagamine 2002). 

Figure 3 shows Color-Magnitude Diagrams 
(CMD) for the simulated and observed galaxies. 
Because the simulation cannot predict the reden- 
ning correction, we plot each simulated galaxy as a 
dotted region (which appears almost like a single 
line), spanning a range in redenning corrections 
from E(B - V) = to E(B - V) = 0.3, consistent 
with values found by Steidel et al. (1999), and 
a range of escape fractions from /esc = 0.1 to 
/esc = 1-0. 

One can see that simulated galaxies have sim- 
ilar colors to the observed ones, which suggests 
that the simulation captures at least some of the 
physics, which is going on in these galaxies. Let us 
now imagine a "super-simulation" , with so large a 
box that it can be considered infinite. Because 
the colors of our brightest galaxies do not change 
much with magnitude, it is plausible to assume 
that simulated galaxies in the magnitude range 
23 < / < 25 would have similar colors to the 
observed ones. If we plausibly assume that the 
simulated luminosity function, when extended to 
higher luminosities, would match the observed one 
reasonably well, then the total luminosity density 
in such a "super-simulation" would be a factor of 
3 higher than in our simulation, and, thus, the 
rcdshift of reionization would be even higher than 
z = 7 - the rcdshift of reionization in our sim- 
ulation - by about Az ~ 1.5. Such a "super- 
simulation" would serve as a strong counter ex- 
ample to claims that LBGs at z ~ 4 are too young 
to reionize the universe by z ~ 6. 

Why, then, is observational determination of 
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Fig. 3. — Color-Magnitude Diagrams for Simu- 
lated Galaxies at z = 4.0 shown as (a) G — I 
bv I and (b) 1Z — I vs I. Each galaxy is repre- 
sented by a region spanning the range of of red- 
dening from E(B - V) = to E(B - V) = 0.3 and 
the range of escape fractions from /esc = 0.1 to 
/esc = 1-0. Squares are observed LBGs from Stei- 
del et al. (1999) in the redshift range from z = 3.5 
to z = 4.8. A cross within the square indicates 
that only a minimum G magnitude could be mea- 
sured. 
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Fig. 4. — Three illustrative star formation histo- 
ries for simulated galaxies with different luminosi- 
ties: / = 26.9 (a), / = 27.2 (6), / = 28.0 (c). 
Despite the very different histories, all three have 
very similar G — I and 1Z — I colors. The tempo- 
ral sampling of star formation histories is 10 times 
higher after 1.2 Gyr, showing fine structure vari- 
ability on a few Myr time scale due to formation 
of individual star clusters, which are resolved (by 
mass) in the simulation. 
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Fig. 5. — Total star formation rates in galaxies 
that appear within a given magnitude bin at z = 
4 as a function of cosmic time. Each indicated 
magnitude is the central value of a bin of width 

1.0 magnitudes. 

ages so misleading? Figure 4 illustrates the diver- 
sity of star formation histories in simulated galax- 
ies. While our brightest galaxies (I < 27, such as 
the one shown in panel a) arc experiencing a peak 
of star formation at around z = 4, one magnitude 
dimmer galaxies (I > 28) were more luminous in 
the past - indeed, at the time of reionization (at 
z ~ 7 in this simulation) - and by z ~ 4 their star 
formation rate has decreased by a factor of sev- 
eral. In fact, the galaxy shown in panel (c) has 

1.1 x 1O 9 M in stars, whereas the galaxy in panel 
(a), while being more than a magnitude brighter, 
has a stellar mass of only 7 x 10 8 M Q . We, there- 
fore, conclude that galaxies that reionized the uni- 
verse at z ~ 7 appear at z ~ 4 within the magni- 
tude range of 28 < I < 30, while LBGs that are 
observed at z ~ 4 with 23 < / < 25 are indeed 
very young, experiencing their first major episode 
of star formation (we refrain from using the word 
"burst" , because the star formation history of a 
galaxy shown in panel a, while strongly rising by 
z ~ 4, can still be hardly called a "burst"). 

Figure 5 illustrates this point further. It shows 
total star formation rates for galaxies that appear 
within different magnitude bins at z = 4. In agree- 
ment with individual examples from Fig. 4, the 
star formation rate at z = 4 is dominated by our 
brightest galaxies, but these galaxies contribute 
only about 20% to the star formation at z ~ 7. 
At that rcdshift more than half of the total star 



Fig. 6. — Mean star formation time for all simu- 
lated galaxies as a function of their magnitudes at 
z = 4. 

formation rate is contributed by galaxies that ap- 
pear within the 28 < I < 30 magnitude bin at 
z ~ 4. 

The large variation in possible star formation 
histories is further confirmed by Figure 6, which 
shows mean mass- weighted star formation ages for 
all simulated galaxies as a function of their mag- 
nitudes at z = 4, 

M * Jo 

where M* is the total stellar mass. Even for 
our brightest galaxies the variation is significant, 
reaching up to a factor of 3 for the dimmest ones. 
A remarkable feature of Fig. 6 is a large concen- 
tration of galaxies with 40 < I < 45 and little 
star formation after 1 Gyr. It is tempting to iden- 
tify these galaxies with dwarf spheroidals: they 
are old and would have luminosities of the order of 
10 5 ~ 6 -Lq at z = 0. A small fraction of these galax- 
ies still have active star formation at z ~ 4, which 
may explain the diverse star formation histories 
of dwarf spheroidals in the Local Group (Grebel 
1998; Mateo 1998). We leave this as a possible 
speculation, as this is not the goal of this paper, 
however. 

Observations at longer wavelengths and/or at 
fainter magnitudes would be more sensitive to 
older populations of stars. Unfortunately, these 
observations are very difficult to make at a red- 
shift of 4, and they will most likely have to wait 
for JWST and 30-meter class ground telescopes. 



7 





-20 -15 -10 -5 
B Absolute Magnitude 



-20 -15 -10 -5 
B Absolute Magnitude 



Fig. 7. — Mass-to-light ratio for simulated galax- 
ies at z = 4. Squares show the luminosity function 
of the simulated galaxies as would be observed in 
the K band with E(B-V) = 0.15 reddening. Ver- 
tical lines through the symbols represent 1-sigma 
Poisson errors. Crosses show the stellar mass-to- 
light ratio in solar units computed with the same 
reddening. The vertical lines on the left show the 
stellar mass-to-light ranges found by Dickinson et 
al. (2003) assuming metallicities of 1/5 solar (dot- 
ted) and solar (dashed). 

But those instruments are on the horizon, and one 
might hope to indeed observe the culprits of reion- 
ization within the 10-year time frame. 

Whatever the global history of star formation 
we would at least expect that the total stellar con- 
tent of the universe would be monotonic with time. 
Thus, we have compared the stellar content of the 
simulated galaxies with the global history of total 
stellar mass as presented in Dickinson et al. (2003). 
In Figure 7 we show the B band rest frame lumi- 
nosity function and mass-to- light ratios of the sim- 
ulated galaxies, calculated from the K band mag- 
nitudes. Vertical lines give the comparison with 
Dickinson et al. (2003) estimates. Our results are 
in comfortable agreement with these estimates for 
all galaxies brighter than dwarf spheroidal candi- 
dates. 

Finally, Figure 8 shows cumulative mass and 
luminosity densities of simulated galaxies as com- 
pared with Dickinson et al. (2003) estimates. 
Again it appears that the simulated galaxies are 
in a reasonable agreement with the observational 
data. We would like to emphasize that we are 
quite comfortable with a factor of 2 agreement: 



Fig. 8. — Cumulative stellar mass and B band lu- 
minosity. Squares show the luminosity function 
of the simulated galaxies at z — 4 as would be 
observed in the K band. Vertical lines represent 
1-sigma Poisson errors. Triangles show the total 
stellar mass in simulated galaxies integrated from 
the lowest magnitude to the plotted magnitude. 
Diamonds trace the luminosity density of simu- 
lated galaxies as would be observed in the K band 
with E(B -V) = 0.15 reddening. The dotted line 
is the mass density range observed by Dickinson 
et al. (2003). The dashed line is their observed 
luminosity density range obtained by integrating 
over a Schechter fit. 

one should remember that the simulation can only 
account for about 1/3 of the total light because 
of the small box size. In fact, in the imaginary 
"super-simulation", discussed above, the total 
stellar mass would be up to a factor of 3 higher (de- 
pending on whether the mass-to-light ratio goes 
down for brighter galaxies) and would not agree 
with Dickinson ct al. (2003). One could hypothe- 
size about a source of this potential discrepancy, 
including anything from a non-standard IMF to 
incorrect stellar synthesis models, but for the pur- 
pose of this paper we are content with the level of 
agreement we obtain, given substantial limitations 
of our simulation. 

It would also be interesting to compare the lu- 
minosities of Lyman Break Galaxies and their dark 
matter masses, because the relation of LBGs and 
their dark matter halos is still poorly understood 
(Somcrville, Primack, & Faber 2001; Wechsler et 
al. 2001; Somcrville 2002; Benson ct al. 2002a,b). 
In Figure 9 we show this relation. As one can see, 
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Fig. 9. — Scatter plot of absolute magnitude at 
1600 A versus total (dark matter + gas + stellar) 
mass of simulated galaxies at redshift 4. Redden- 
ing equivalent to a color excess of 0.15 has been 
applied as well as extinction due to absorption by 
the Lyman alpha forest. 

Lyman Break Galaxies represent a diverse popu- 
lation. For example, our brightest galaxies with 
absolute magnitudes between -18 and -17 cover 
about two magnitudes in the total mass: some of 
them are the most massive galaxies in the simu- 
lations, but some of them have significantly lower 
masses and are experiencing a starburst at z = 4. 
There exists a general trend of less massive galax- 
ies being dimmer, but the scatter around the mean 
relation is very large. 

4. Discussion 

So, how well does the simulation do? While the 
computational box of our simulation is clearly in- 
sufficient to have even a single galaxy as luminous 
as those of Steidel et al. (1999), simulated galaxies 
do appear as fainter cousins of the observed ones. 
They have similar colors, and simulated and ob- 
served galaxies together can be fitted by a single 
luminosity function. 

They are also faring well in their abilities to 
reionize the universe. Star formation histories of 
simulated galaxies are diverse, and vary systemat- 
ically with magnitude. The brightest galaxies at 
z = 4 are indeed quite young, in accord with con- 
clusions of Shapely et al. (2001) and Ferguson et 
al. (2002). But these brightest galaxies contribute 
only about 25% of ionizing photons at z > 6. 



The bulk of work of reionizing the universe was 
done by dimmer galaxies, those that fall within 
the 28 < I < 30 magnitude range at z ~ 4, but 
which arc not necessarily less massive than the 
bright ones. 

We expect that detailed simulations will be- 
come increasingly important in sorting out the 
various possibilities. The general picture that our 
simulation produces is at least consistent with ob- 
servations of Lyman Break Galaxies, an encourag- 
ing message to theorists. 

The recent results from the WMAP mission 
(Kogut et al. 2003) indicate that the optical depth 
to Thompson scattering is significantly larger than 
in our simulation. While this measurement cannot 
pin down the redshift of reionization, it suggests, 
if correct, that there was a considerable ionizing 
flux at early times. Our simulation provides no 
support for a stellar origin of this ionizing flux. 
Perhaps, very massive population III stars, which 
we have not considered, could be responsible. 

Because, as we have discussed, the photoion- 
ization feedback likely plays a significant role in 
the evolution of faint Lyman Break Galaxies, one 
might wonder how the WMAP results affect the 
conclusions of this paper. As has been shown by 
several previous studies (Gnedin 2000b and ref- 
erences herein), the photoionization feedback re- 
sults in the substantial loss of the gas mass of ob- 
jects with circular velocities below about 40km/s. 
However, this characteristic circular velocity is es- 
sentially redshift independent. So even if the uni- 
verse was fully ionized throughout its entire his- 
tory, the effect of the photoionization feedback at 
z = 4 would be essentially at the same level as in 
our simulations. Our conclusion about a reason- 
able agreement between the simulation and the 
data is, in fact, independent of the redshift of 
reionization as long as it is above about 6. 

We are grateful to Chuck Steidel for provid- 
ing us with his custom filter shapes and for en- 
lightening comments during his visit to CU. This 
work was supported by NSF grant AST-0134373. 
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tional Computational Science Alliance under grant 
AST-960015N and utilized the SGI/CRAY Origin 
2000 array at the National Center for Supercom- 
puting Applications (NCSA). 
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